R Notebook

The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file

Setup

Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd

  |                                                                                                                                                                                                                                                                                             
  |                                                                                                                                                                                                                                                                                       |   0%
  |                                                                                                                                                                                                                                                                                             
  |.............                                                                                                                                                                                                                                                                          |   5%
  |                                                                                                                                                                                                                                                                                             
  |...........................                                                                                                                                                                                                                                                            |  10%
  |                                                                                                                                                                                                                                                                                             
  |........................................                                                                                                                                                                                                                                               |  14%
  |                                                                                                                                                                                                                                                                                             
  |.....................................................                                                                                                                                                                                                                                  |  19%
  |                                                                                                                                                                                                                                                                                             
  |..................................................................                                                                                                                                                                                                                     |  24%
  |                                                                                                                                                                                                                                                                                             
  |................................................................................                                                                                                                                                                                                       |  29%
  |                                                                                                                                                                                                                                                                                             
  |.............................................................................................                                                                                                                                                                                          |  33%
  |                                                                                                                                                                                                                                                                                             
  |..........................................................................................................                                                                                                                                                                             |  38%
  |                                                                                                                                                                                                                                                                                             
  |........................................................................................................................                                                                                                                                                               |  43%
  |                                                                                                                                                                                                                                                                                             
  |.....................................................................................................................................                                                                                                                                                  |  48%
  |                                                                                                                                                                                                                                                                                             
  |..................................................................................................................................................                                                                                                                                     |  52%
  |                                                                                                                                                                                                                                                                                             
  |...............................................................................................................................................................                                                                                                                        |  57%
  |                                                                                                                                                                                                                                                                                             
  |.............................................................................................................................................................................                                                                                                          |  62%
  |                                                                                                                                                                                                                                                                                             
  |..........................................................................................................................................................................................                                                                                             |  67%
  |                                                                                                                                                                                                                                                                                             
  |.......................................................................................................................................................................................................                                                                                |  71%
  |                                                                                                                                                                                                                                                                                             
  |.....................................................................................................................................................................................................................                                                                  |  76%
  |                                                                                                                                                                                                                                                                                             
  |..................................................................................................................................................................................................................................                                                     |  81%
  |                                                                                                                                                                                                                                                                                             
  |...............................................................................................................................................................................................................................................                                        |  86%
  |                                                                                                                                                                                                                                                                                             
  |............................................................................................................................................................................................................................................................                           |  90%
  |                                                                                                                                                                                                                                                                                             
  |..........................................................................................................................................................................................................................................................................             |  95%
  |                                                                                                                                                                                                                                                                                             
  |.......................................................................................................................................................................................................................................................................................| 100%

output file: /tmp/RtmpPXtcTe/file315a4f0aaec

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose


Attaching package: ‘ftExtra’

The following object is masked from ‘package:flextable’:

    separate_header

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

nonSpikes %>% head

raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180   3-1-S-20    3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20    3-2-B-5  3-2-B-500   3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53 
  4.382762   5.084370   4.678155   5.926793   5.237899   3.867372   5.596986   4.563925   4.967078   4.765325   5.298189   4.815413   4.407171   4.665767   4.535506   5.052774   5.446895   5.077076   4.960559   3.757001   4.947850   4.690475   4.971215   5.197389   4.867073   4.144527 
 3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500   3-3-B-53  3-3-S-180   3-3-S-20  3-3-S-500   3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500   4-3-O-53  4-3-S-0-2  4-3-S-180   4-3-S-20 
  4.317769   4.955861   3.441315   5.775953   5.283504   5.340186   5.501008   5.042675   5.026191   4.920511   4.302272   4.376545   4.946276   4.475293   4.367380   4.652247   4.290391   4.321040   5.034992   4.637150   5.244157   4.021406   4.520193   2.805419   4.545389   4.648545 
 4-3-S-500   4-3-S-53  5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53  5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500   5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2  C_5P1B_20 C_5P1B_500  C_5P1B_53 
  4.747469   4.443062   4.354328   4.570223   4.145389   4.095409   4.085605   4.167853   4.681186   5.103463   5.441884   5.083902   4.899043   4.294475   4.887705   4.875368   4.208047   4.295552   4.779480   4.864873   5.423950   4.773348   5.107284 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2   3-1-S-180    3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2   3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500    3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5 
0.011295779 0.008105812 0.012287837 0.003007367 0.005702428 0.042036654 0.008402445 0.010684687 0.007005046 0.008998109 0.006977728 0.009101954 0.014352266 0.009499526 0.014455797 0.006241515 0.007677090 0.008637849 0.012140483 0.017485268 0.008113446 0.010226522 0.006023986 0.009093441 
  3-2-S-500    3-2-S-53   3-3-B-0-2   3-3-B-1-2   3-3-B-180    3-3-B-20     3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500    3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5   4-3-B-500    4-3-B-53   4-3-O-1-2   4-3-O-180     4-3-O-5   4-3-O-500 
0.009069482 0.014334723 0.012184340 0.009141746 0.064686366 0.004660697 0.007687693 0.005758500 0.006703955 0.008822806 0.008114157 0.009782327 0.012208135 0.010270972 0.008396074 0.009961698 0.009821496 0.008589214 0.008862728 0.011374560 0.008911489 0.009510012 0.007729699 0.011732768 
   4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53   5-1-S-1-2   5-1-S-180    5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53   5-5-B-0-2   5-5-B-180     5-5-B-5   5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2  C_5P1B_180  C_5P1B_1P2 
0.009779128 0.124685269 0.012987829 0.009992095 0.013958218 0.014609567 0.010345631 0.012334306 0.014569903 0.017536336 0.012430848 0.015666936 0.007652523 0.008220342 0.005410565 0.009486091 0.009221795 0.016242008 0.011576142 0.009351730 0.015098839 0.009240293 0.007351504 0.005571680 
  C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.003605969 0.007392898 0.005904266 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied observed species numbers

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-224.389  -40.954    2.457   53.730  199.561 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          292239.870 138034.341   2.117 0.037854 *  
log(Size_Class)          32.107      8.363   3.839 0.000271 ***
I(log(Size_Class)^2)     -6.364      1.556  -4.090 0.000115 ***
lat                  -15220.663   7192.791  -2.116 0.037947 *  
I(lat^2)                198.241     93.636   2.117 0.037856 *  
depth                     5.152      3.312   1.555 0.124459    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 78.62 on 69 degrees of freedom
Multiple R-squared:  0.2554,    Adjusted R-squared:  0.2015 
F-statistic: 4.734 on 5 and 69 DF,  p-value: 0.0008913

Richness

Rarified chao1 estimates

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-539.32 -131.99  -38.59  128.05 1210.56 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)           1.071e+06  4.831e+05   2.217  0.02989 * 
log(Size_Class)       8.568e+01  2.927e+01   2.927  0.00463 **
I(log(Size_Class)^2) -1.801e+01  5.446e+00  -3.307  0.00150 **
lat                  -5.582e+04  2.517e+04  -2.218  0.02988 * 
I(lat^2)              7.272e+02  3.277e+02   2.219  0.02978 * 
depth                 2.227e+01  1.159e+01   1.921  0.05884 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 275.2 on 69 degrees of freedom
Multiple R-squared:  0.1982,    Adjusted R-squared:  0.1401 
F-statistic: 3.411 on 5 and 69 DF,  p-value: 0.008206

As above but without latitude and depth

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-466.78 -164.77  -27.36  117.52 1308.99 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           568.437     50.075  11.352  < 2e-16 ***
log(Size_Class)        86.353     29.738   2.904  0.00489 ** 
I(log(Size_Class)^2)  -18.426      5.532  -3.331  0.00137 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 279.9 on 72 degrees of freedom
Multiple R-squared:  0.1341,    Adjusted R-squared:  0.1101 
F-statistic: 5.576 on 2 and 72 DF,  p-value: 0.005603

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.35466 -0.22941  0.03062  0.30917  0.76763 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.683e+03  8.051e+02   2.090  0.04028 *  
log(Size_Class)       1.990e-01  4.878e-02   4.079  0.00012 ***
I(log(Size_Class)^2) -3.661e-02  9.075e-03  -4.034  0.00014 ***
lat                  -8.744e+01  4.195e+01  -2.084  0.04083 *  
I(lat^2)              1.139e+00  5.461e-01   2.085  0.04079 *  
depth                 2.099e-02  1.932e-02   1.087  0.28095    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4586 on 69 degrees of freedom
Multiple R-squared:  0.2863,    Adjusted R-squared:  0.2346 
F-statistic: 5.537 on 5 and 69 DF,  p-value: 0.000241

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.014482 -0.005102 -0.002446  0.000812  0.099932 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.506e+01  2.657e+01  -0.943   0.3490  
log(Size_Class)      -3.808e-03  1.610e-03  -2.365   0.0208 *
I(log(Size_Class)^2)  6.620e-04  2.996e-04   2.210   0.0304 *
lat                   1.305e+00  1.385e+00   0.942   0.3493  
I(lat^2)             -1.697e-02  1.803e-02  -0.941   0.3498  
depth                -3.651e-04  6.377e-04  -0.573   0.5688  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01514 on 69 degrees of freedom
Multiple R-squared:  0.09266,   Adjusted R-squared:  0.02691 
F-statistic: 1.409 on 5 and 69 DF,  p-value: 0.2318

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20       21       22       23       24       25       26       27       28       29       30       31       32 
231.2129 231.2129 197.7419 197.7419 213.1276 161.0777 161.0777 193.7155 203.6695 305.0142 305.0142 271.5432 271.5432 286.9290 234.8790 234.8790 267.5169 267.5169 334.5602 334.5602 301.0891 301.0891 316.4749 264.4250 264.4250 297.0628 307.0168 307.0168 338.4386 338.4386 304.9676 304.9676 
      33       34       35       36       37       38       39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60       61       62       63       64 
320.3533 320.3533 268.3034 268.3034 300.9413 300.9413 326.5221 293.0510 293.0510 308.4368 308.4368 256.3869 256.3869 256.3869 289.0247 289.0247 298.9787 298.9787 294.4742 294.4742 261.0032 261.0032 276.3890 276.3890 224.3390 224.3390 224.3390 256.9769 256.9769 266.9308 266.9308 253.1023 
      65       66       67       68       69       70       71       72       73       74       75 
219.6313 219.6313 235.0171 235.0171 182.9671 182.9671 182.9671 215.6050 215.6050 225.5589 225.5589 

$se.fit
 [1] 27.50522 27.50522 26.75531 26.75531 26.61253 28.31263 28.31263 29.90239 34.68922 19.29454 19.29454 18.66638 18.66638 17.52990 20.42368 20.42368 21.91454 21.91454 19.62244 19.62244 19.15178 19.15178 17.56592 20.47128 20.47128 21.61090 28.71618 28.71618 19.62022 19.62022 19.11805
[32] 19.11805 17.28640 17.28640 19.98620 19.98620 20.96913 20.96913 18.91721 18.26594 18.26594 16.30436 16.30436 18.79606 18.79606 18.79606 19.80604 19.80604 27.06275 27.06275 19.68177 19.68177 18.77782 18.77782 16.99969 16.99969 18.77935 18.77935 18.77935 19.87208 19.87208 26.61473
[63] 26.61473 24.92790 23.95655 23.95655 22.77274 22.77274 23.58385 23.58385 23.58385 24.60241 24.60241 29.78621 29.78621

$df
[1] 69

$residual.scale
[1] 78.62265

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20       21       22       23       24       25       26       27       28       29       30       31       32 
480.6067 480.6067 353.4462 353.4462 448.7195 264.4622 264.4622 406.6256 382.2290 680.1825 680.1825 553.0219 553.0219 648.2953 464.0380 464.0380 606.2014 606.2014 756.4058 756.4058 629.2453 629.2453 724.5186 540.2614 540.2614 682.4247 658.0281 658.0281 760.2027 760.2027 633.0422 633.0422 
      33       34       35       36       37       38       39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60       61       62       63       64 
728.3155 728.3155 544.0583 544.0583 686.2216 686.2216 721.4327 594.2722 594.2722 689.5455 689.5455 505.2883 505.2883 505.2883 647.4517 647.4517 623.0550 623.0550 624.4071 624.4071 497.2465 497.2465 592.5199 592.5199 408.2626 408.2626 408.2626 550.4260 550.4260 526.0294 526.0294 502.0359 
      65       66       67       68       69       70       71       72       73       74       75 
374.8753 374.8753 470.1487 470.1487 285.8914 285.8914 285.8914 428.0548 428.0548 403.6582 403.6582 

$se.fit
 [1]  96.26104  96.26104  93.63656  93.63656  93.13686  99.08674  99.08674 104.65049 121.40314  67.52582  67.52582  65.32740  65.32740  61.35005  71.47751  71.47751  76.69512  76.69512  68.67338  68.67338  67.02618  67.02618  61.47610  71.64408  71.64408  75.63246 100.49907 100.49907
[29]  68.66562  68.66562  66.90815  66.90815  60.49784  60.49784  69.94644  69.94644  73.38642  73.38642  66.20526  63.92598  63.92598  57.06097  57.06097  65.78128  65.78128  65.78128  69.31593  69.31593  94.71251  94.71251  68.88102  68.88102  65.71742  65.71742  59.49445  59.49445
[57]  65.72280  65.72280  65.72280  69.54705  69.54705  93.14456  93.14456  87.24108  83.84161  83.84161  79.69860  79.69860  82.53725  82.53725  82.53725  86.10195  86.10195 104.24388 104.24388

$df
[1] 69

$residual.scale
[1] 275.1586

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19       20       21       22       23       24       25       26       27       28       29       30       31       32 
4.524800 4.524800 4.344339 4.344339 4.334988 4.022197 4.022197 4.154170 4.370926 4.974911 4.974911 4.794451 4.794451 4.785100 4.472309 4.472309 4.604282 4.604282 5.165240 5.165240 4.984780 4.984780 4.975429 4.662638 4.662638 4.794611 5.011367 5.011367 5.207335 5.207335 5.026874 5.026874 
      33       34       35       36       37       38       39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57       58       59       60       61       62       63       64 
5.017523 5.017523 4.704732 4.704732 4.836705 4.836705 5.152692 4.972232 4.972232 4.962881 4.962881 4.650090 4.650090 4.650090 4.782063 4.782063 4.998819 4.998819 4.985785 4.985785 4.805324 4.805324 4.795974 4.795974 4.483183 4.483183 4.483183 4.615156 4.615156 4.831912 4.831912 4.762373 
      65       66       67       68       69       70       71       72       73       74       75 
4.581912 4.581912 4.572561 4.572561 4.259770 4.259770 4.259770 4.391743 4.391743 4.608499 4.608499 

$se.fit
 [1] 0.1604232 0.1604232 0.1560494 0.1560494 0.1552166 0.1651324 0.1651324 0.1744046 0.2023236 0.1125347 0.1125347 0.1088710 0.1088710 0.1022425 0.1191204 0.1191204 0.1278158 0.1278158 0.1144472 0.1144472 0.1117021 0.1117021 0.1024526 0.1193980 0.1193980 0.1260448 0.1674861 0.1674861
[29] 0.1144342 0.1144342 0.1115054 0.1115054 0.1008223 0.1008223 0.1165688 0.1165688 0.1223017 0.1223017 0.1103340 0.1065354 0.1065354 0.0950946 0.0950946 0.1096274 0.1096274 0.1096274 0.1155180 0.1155180 0.1578425 0.1578425 0.1147932 0.1147932 0.1095210 0.1095210 0.0991501 0.0991501
[57] 0.1095299 0.1095299 0.1095299 0.1159032 0.1159032 0.1552295 0.1552295 0.1453911 0.1397257 0.1397257 0.1328212 0.1328212 0.1375519 0.1375519 0.1375519 0.1434927 0.1434927 0.1737270 0.1737270

$df
[1] 69

$residual.scale
[1] 0.4585638

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9          10          11          12          13          14          15          16          17          18          19          20          21          22          23          24 
0.019119982 0.019119982 0.021726336 0.021726336 0.020969819 0.024752823 0.024752823 0.022023135 0.018632849 0.010604187 0.010604187 0.013210541 0.013210541 0.012454023 0.016237028 0.016237028 0.013507339 0.013507339 0.006862322 0.006862322 0.009468676 0.009468676 0.008712158 0.012495163 
         25          26          27          28          29          30          31          32          33          34          35          36          37          38          39          40          41          42          43          44          45          46          47          48 
0.012495163 0.009765474 0.006375189 0.006375189 0.005809293 0.005809293 0.008415647 0.008415647 0.007659129 0.007659129 0.011442134 0.011442134 0.008712445 0.008712445 0.006592059 0.009198414 0.009198414 0.008441896 0.008441896 0.012224901 0.012224901 0.012224901 0.009495212 0.009495212 
         49          50          51          52          53          54          55          56          57          58          59          60          61          62          63          64          65          66          67          68          69          70          71          72 
0.006104927 0.006104927 0.009352450 0.009352450 0.011958804 0.011958804 0.011202287 0.011202287 0.014985291 0.014985291 0.014985291 0.012255603 0.012255603 0.008865318 0.008865318 0.013176838 0.015783192 0.015783192 0.015026675 0.015026675 0.018809679 0.018809679 0.018809679 0.016079991 
         73          74          75 
0.016079991 0.012689706 0.012689706 

$se.fit
 [1] 0.005295216 0.005295216 0.005150846 0.005150846 0.005123358 0.005450654 0.005450654 0.005756710 0.006678255 0.003714522 0.003714522 0.003593590 0.003593590 0.003374800 0.003931900 0.003931900 0.004218915 0.004218915 0.003777648 0.003777648 0.003687038 0.003687038 0.003381734
[24] 0.003941063 0.003941063 0.004160460 0.005528345 0.005528345 0.003777221 0.003777221 0.003680545 0.003680545 0.003327921 0.003327921 0.003847678 0.003847678 0.004036908 0.004036908 0.003641880 0.003516499 0.003516499 0.003138862 0.003138862 0.003618557 0.003618557 0.003618557
[47] 0.003812994 0.003812994 0.005210033 0.005210033 0.003789070 0.003789070 0.003615044 0.003615044 0.003272725 0.003272725 0.003615340 0.003615340 0.003615340 0.003825708 0.003825708 0.005123781 0.005123781 0.004799037 0.004612036 0.004612036 0.004384134 0.004384134 0.004540285
[70] 0.004540285 0.004540285 0.004736375 0.004736375 0.005734343 0.005734343

$df
[1] 69

$residual.scale
[1] 0.01513617

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
  bold(i = ~ p< 0.05, j = "p") %>%
  colformat_md() %>%
  set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Metric	Term	Estimate	Standard Error	T Value	p
Observed ASVs	Intercept	2.9 x 105	1.4 x 105	2.12	0.038
	log(Size Class)	3.2 x 101	8.4 x 100	3.84	< 0.001
	log(Size Class)2	-6.4 x 100	1.6 x 100	-4.09	< 0.001
	Latitude	-1.5 x 104	7.2 x 103	-2.12	0.038
	Latitude2	2.0 x 102	9.4 x 101	2.12	0.038
	Depth	5.2 x 100	3.3 x 100	1.56	0.124
Richness (Chao1)	Intercept	1.1 x 106	4.8 x 105	2.22	0.030
	log(Size Class)	8.6 x 101	2.9 x 101	2.93	0.005
	log(Size Class)2	-1.8 x 101	5.4 x 100	-3.31	0.001
	Latitude	-5.6 x 104	2.5 x 104	-2.22	0.030
	Latitude2	7.3 x 102	3.3 x 102	2.22	0.030
	Depth	2.2 x 101	1.2 x 101	1.92	0.059
Diversity (Shannon H)	Intercept	1.7 x 103	8.1 x 102	2.09	0.040
	log(Size Class)	2.0 x 10-1	4.9 x 10-2	4.08	< 0.001
	log(Size Class)2	-3.7 x 10-2	9.1 x 10-3	-4.03	< 0.001
	Latitude	-8.7 x 101	4.2 x 101	-2.08	0.041
	Latitude2	1.1 x 100	5.5 x 10-1	2.08	0.041
	Depth	2.1 x 10-2	1.9 x 10-2	1.09	0.281
Evenness (Pielou J)	Intercept	-2.5 x 101	2.7 x 101	-0.94	0.349
	log(Size Class)	-3.8 x 10-3	1.6 x 10-3	-2.37	0.021
	log(Size Class)2	6.6 x 10-4	3.0 x 10-4	2.21	0.030
	Latitude	1.3 x 100	1.4 x 100	0.94	0.349
	Latitude2	-1.7 x 10-2	1.8 x 10-2	-0.94	0.350
	Depth	-3.7 x 10-4	6.4 x 10-4	-0.57	0.569

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013926 -0.005094 -0.002507  0.000907  0.105946 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.955e+01  2.753e+01  -0.710   0.4800  
log(Size_Class)      -3.291e-03  1.668e-03  -1.974   0.0524 .
I(log(Size_Class)^2)  5.747e-04  3.103e-04   1.852   0.0683 .
lat                   1.017e+00  1.434e+00   0.709   0.4809  
I(lat^2)             -1.321e-02  1.867e-02  -0.707   0.4818  
depth                -2.353e-04  6.605e-04  -0.356   0.7228  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01568 on 69 degrees of freedom
Multiple R-squared:  0.06733,   Adjusted R-squared:  -0.0002567 
F-statistic: 0.9962 on 5 and 69 DF,  p-value: 0.4267

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7             8             9            10            11            12            13            14            15            16            17            18            19            20 
 0.0116882006  0.0116882006  0.0135742540  0.0135742540  0.0133932226  0.0160282544  0.0160282544  0.0140111437  0.0099141263  0.0043215048  0.0043215048  0.0062075582  0.0062075582  0.0060265267  0.0086615586  0.0086615586  0.0066444478  0.0066444478  0.0010939134  0.0010939134 
           21            22            23            24            25            26            27            28            29            30            31            32            33            34            35            36            37            38            39            40 
 0.0029799668  0.0029799668  0.0027989353  0.0054339672  0.0054339672  0.0034168564 -0.0006801609 -0.0006801609  0.0002000096  0.0002000096  0.0020860630  0.0020860630  0.0019050316  0.0019050316  0.0045400634  0.0045400634  0.0025229527  0.0025229527  0.0008938114  0.0027798648 
           41            42            43            44            45            46            47            48            49            50            51            52            53            54            55            56            57            58            59            60 
 0.0027798648  0.0025988333  0.0025988333  0.0052338652  0.0052338652  0.0052338652  0.0032167544  0.0032167544 -0.0008802629 -0.0008802629  0.0033080946  0.0033080946  0.0051941480  0.0051941480  0.0050131165  0.0050131165  0.0076481484  0.0076481484  0.0076481484  0.0056310377 
           61            62            63            64            65            66            67            68            69            70            71            72            73            74            75 
 0.0056310377  0.0015340203  0.0015340203  0.0066431363  0.0085291897  0.0085291897  0.0083481582  0.0083481582  0.0109831901  0.0109831901  0.0109831901  0.0089660794  0.0089660794  0.0048690620  0.0048690620 

$se.fit
 [1] 0.005484918 0.005484918 0.005335376 0.005335376 0.005306903 0.005645925 0.005645925 0.005962946 0.006917505 0.003847596 0.003847596 0.003722331 0.003722331 0.003495703 0.004072762 0.004072762 0.004370059 0.004370059 0.003912983 0.003912983 0.003819127 0.003819127 0.003502885
[24] 0.004082253 0.004082253 0.004309509 0.005726399 0.005726399 0.003912541 0.003912541 0.003812401 0.003812401 0.003447144 0.003447144 0.003985522 0.003985522 0.004181531 0.004181531 0.003772351 0.003642478 0.003642478 0.003251313 0.003251313 0.003748193 0.003748193 0.003748193
[47] 0.003949595 0.003949595 0.005396683 0.005396683 0.003924814 0.003924814 0.003744554 0.003744554 0.003389971 0.003389971 0.003744860 0.003744860 0.003744860 0.003962765 0.003962765 0.005307342 0.005307342 0.004970964 0.004777263 0.004777263 0.004541196 0.004541196 0.004702942
[70] 0.004702942 0.004702942 0.004906057 0.004906057 0.005939777 0.005939777

$df
[1] 69

$residual.scale
[1] 0.01567843

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

Metric	Term	Estimate	Standard Error	T Value	p
Richness (Breakaway – Betta)	Intercept	7.1 x 106	2.4 x 102	NA	< 0.001
	log(Size Class)	1.2 x 102	6.1 x 101	NA	0.058
	log(Size Class)2	-5.0 x 101	1.2 x 101	NA	< 0.001
	Latitude	-3.7 x 105	6.1 x 100	NA	< 0.001
	Latitude2	4.8 x 103	1.6 x 10-1	NA	< 0.001
	Depth	1.5 x 102	1.0 x 101	NA	< 0.001
Richness (Breakaway – LM)	Intercept	7.1 x 106	3.3 x 106	2.13	0.036
	log(Size Class)	2.4 x 102	2.0 x 102	1.21	0.231
	log(Size Class)2	-7.5 x 101	3.8 x 101	-2.00	0.050
	Latitude	-3.7 x 105	1.7 x 105	-2.13	0.037
	Latitude2	4.8 x 103	2.3 x 103	2.13	0.037
	Depth	1.5 x 102	8.0 x 101	1.89	0.064
Diversity (Shannon H)	Intercept	1.7 x 103	8.1 x 102	2.09	0.040
	log(Size Class)	2.0 x 10-1	4.9 x 10-2	4.08	< 0.001
	log(Size Class)2	-3.7 x 10-2	9.1 x 10-3	-4.03	< 0.001
	Latitude	-8.7 x 101	4.2 x 101	-2.08	0.041
	Latitude2	1.1 x 100	5.5 x 10-1	2.08	0.041
	Depth	2.1 x 10-2	1.9 x 10-2	1.09	0.281
Evenness (Pielou J)	Intercept	-2.0 x 101	2.8 x 101	-0.71	0.480
	log(Size Class)	-3.3 x 10-3	1.7 x 10-3	-1.97	0.052
	log(Size Class)2	5.7 x 10-4	3.1 x 10-4	1.85	0.068
	Latitude	1.0 x 100	1.4 x 100	0.71	0.481
	Latitude2	-1.3 x 10-2	1.9 x 10-2	-0.71	0.482
	Depth	-2.4 x 10-4	6.6 x 10-4	-0.36	0.723

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)